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Smoothed particle hydrodynamics (SPH) discretization techniques are gen- 
eralized to develop a method, smoothed particle interpolation (SPI), for solving 
initial value problems of systems of non-hydrodynamical nature. Under this 
approach, SPH is viewed as strictly an interpolation scheme and, as such, suit- 
able for solving general hyperbolic and parabolic equations. The SPI method is 
tested on (1) the wave equation with inhomogeneous sound speed and (2) Burgers 
equation. The efficiency of SPI is studied by comparing SPI solutions to those 
obtained with standard finite difference methods. It is shown that the power of 
SPI arises when the smoothing particles are free to move. 
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1. INTRODUCTION 



Perhaps one of the most powerful aspects in computational studies is the development of 
a fully adaptive scheme for solving partial differential equations. Without adaptive numeri- 
cal schemes, current large memory and fast supercomputers are not capable yet of obtaining 
highly accurate solutions to physically interesting three-dimensional problems. For this rea- 
son, an optimal management of computational resources has been one of the driving forces 
in the effort to achieve both temporal and spatial adaptability as a way of obtaining a higher 
degree of accuracy. 

The goal for reaching optimal temporal adaptability consists of designing algorithms 
that evolve each region of the computational domain with time-steps according to their 
local time-scales. In particle methods, local time-step algorithms are aimed at subsets of 
particles. The main difficulty here arises in designing a time evolution algorithm compatible 
with the neighbor-finding scheme (tree-code). In the case of finite differencing, the target of 
adaptability are subgrids. Each subgrid, for instance, is evolved according to its own Courant 
condition. Special attention must be paid at the subgrid boundaries where gridpoints in 
coarser meshes are not updated as often as in the finer meshes. Other approaches to temporal 
adaptability exploit multiple time-stepping applied at the equation instead of the domain 
level. These methods rely on the possibility of identifying in the evolution equations the 
slowly varying terms that require time-step At in contrast to rapidly varying terms with 
time-step St < At. The slowly varying terms are then computed once every large time-step 
At, and the rapidly varying terms are subcycled At/5t times for each large time-step At. 
An example of this approach is called acoustic-subcycling, which has successfully improved 
the efficiency of compressible flow codes (O'Rourke & Amsden 1986) when applied to low 
Mach number problems. 

Technically speaking, spatial adaptability exhibits a higher degree of difficulty than tem- 
poral one. For finite difference methods, this translates into designing an efficient remeshing 
algorithm. Spatially adaptive finite difference methods for hyperbolic equations have been 
proposed by Berger and Oliger (1984) using multiple component grids. This approach was 
the keystone for Choptuik's (1993) finding of scaling and critical behavior in the gravita- 
tional collapse of a scalar field. Adaptive finite element methods, which have been mainly 
used for elliptic problems, have been also applied to hyperbolic equations (Davis & Flaherty 
1982). Finally, particle methods have been explicitly designed to exploit adaptability; a set 
of moving points carry the information of the dynamics of the system. Smoothed particle 
hydrodynamics (SPH) (for a review, see Benz 1990) is an example of such methods and 
constitutes the focus of this paper. 

SPH has become a remarkably robust hydrodynamic method, in particular for astrophys- 
ical systems where vast range of scales are present. It has been successfully applied to stellar 
encounters(Davies, Benz, & Hills 1991), i galaxy formation (Hernquist & Katz 1989), tidal 
disruptions (Evans & Kochanek 1989; Laguna, Miller, Zurek, & Davies 1993), and more re- 
cently to problems outside the astrophysical domain, such as hypervelocity impacts. SPH has 
also been reformulated to describe relativistic fluids in curved space-times (Laguna, Miller 
& Zurek 1993). SPH has been also used in combination with other methods in studying 
mixed systems where additional physics is required (hydrodynamics + magnetism, hydro- 
dynamics + gravity, hydrodynamics + collisionless matter, etc.). In most of these cases, 
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SPH is only applied to the hydro dynamic equations, and other numerical methods are used 
in the remaining set of equations. Mann (1993) has solved the relativistic, self-gravitating 
spherical collapse using finite element methods to compute the gravitational potentials and 
SPH for the hydrodynamics equations. In the study of large scale structure formation, SPH 
is combined with iV-body methods to simulate hybrid systems of Hot (SPH) + Cold (N- 
body) Dark Matter. There are, however, examples of studies in which, still in the presence 
of fluids, SPH has been applied to non-hydrodynamical equations. For instance, in the low 
frequency and infinite conductivity limit Stellingwerf and Peterkin (1989) have developed a 
simple and elegant method for smoothed particled magnetohydrodynamics. 

The main goal of this paper is to address the following question: Can the success of 
SPH be extended to purely non-hydrodynamical problems? In other words, can one take 
advantage of SPH gridless nature of computing spatial derivatives and move particles with 
arbitrary "grid" velocity? We present one- dimensional results that point in the direction of 
an affirmative answer to this question. In Section 2 of the paper, we present the method for 
applying SPH discretization techniques to general hyperbolic and parabolic equations, giving 
special attention to truncation errors. Section 3 contains one-dimensional tests consisting of 
solutions to the wave equation with inhomogeneous sound speed; for these tests the particles 
are "fixed" in time (Eulerian framework). In Section 4, we repeat the tests of Section 3, but 
now we allow for particle motion (Lagrangian framework). Solutions to Burgers equation 
using SPH-type methods are obtained in Section 5. Finally, discussion and conclusions are 
given in Section 5. 



The key ingredients of SPH are: (1) the fluid elements of the system are represented 
by particles, (2) the dynamics of the particles is governed by the hydrodynamic equations 
written in the Lagrangian form, and (3) the spatial derivatives at the particle positions are 
calculated via an interpolation kernel. Given a physical function /(x) (pressure, density, 
etc.), its mean smoothed value can be obtained from 



where W(x.,x.';h) is the kernel and h the smoothing length. The kernel W(x, x!;h) has 
compact support in a region of 0(h), and for hydrodynamic calculations the kernel is assumed 
to be of class C 1 , i.e., its first derivatives exist and are continuous. One can show that this 
smoothing procedure is second-order accurate; that is, (/(x)) = /(x) + 0(h 2 ). In particular, 
for the case of an spherically symmetric kernel the truncation errors due to the smoothing 
are 



2. SPI: SMOOTHED PARTICLE INTERPOLATION 




(2.1) 



(f(x)) = f(x) + ah 2 Af + 0(h s ), 



(2.2) 



where 




(2.3) 
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is independent of h; in the one- dimensional case, a = 1/4. The consistency requirement that 
(/(x)) — > /(x) as h — ► 0, imposes the following normalization condition on W: 

[w(x,x';h)d 3 x = l. (2.4) 



The smoothing procedure Q2.1|) provides a natural prescription for calculating spatial deriva- 
tives at the particle locations. One has that, when ignoring surface terms, integration by 
parts of 



< V/(x)> = J W(x, x'; h)Vf{x!)d 3 x' (2.5) 
yields, 

<V/(x)> = / /(x') W(x, x'; /i)dV, (2.6) 



where V and V are gradients with respect to x' and x, respectively. Equation (|2.6f) is 
just the statement that (V/(x)) = V(/(x)). Similarly, it is not difficult to show that for a 
vector field v(x), the smoothed value of its divergence is given by (V • v(x)) = V • (v(x)), 
or equivalently by 

(V • v(x)) = J v(x') • Vlf(x, x'; h)d 3 x'. (2.7) 

Finally, the smoothing of V 2 / is obtained by successively applying properties (V ■ v(x)) = 
V • (v(x)) and (V/(x)) = V(/(x)). One gets that (V 2 /) = V 2 (/), which has as integral 
representation 

(V 2 /(x Q )> = J /(x') V 2 ^(x, x'; h)d 3 x'. (2.8) 

One can show that the second-order accurate nature of the smoothing procedure is pre- 
served under differentiation. The differentiation rules (SH) and ( |2.7| ) to compute (V/(x)) 
and (V- v(x)), together with definition (|2TT|), is all that is needed to derive the SPH equations 
in their integral form (Benz 1990). 

The next step is to obtain the discrete version of Eqs. ( |2.1| ) and (|2.6|)-( |2~8"D ; that is, a 
discrete representation of the integrals in those equations. In order to accomplish this, one 
defines first a number density distribution 

N 

n(x) = £5(x-x ), (2.9) 

a=l 

with {xo}a=i,..,jv the collection of iV-points (particles) where the physical functions are 
known. The integrals in Eqs. ( |2.1| ) and ( |2 . 6| ) - ( PT5| ) are evaluated by multiplying the inte- 
grands by n(x)/(n(x)) = 1 + 0(h 2 ), yielding 

N fix*.) 

[/(x a )] = £ /L^^( Xa , x b ; h) , (2.10) 
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[V/(x )] ^/nrV^fx.x^), (2.11) 
[V ■ v(x a )] = £ • V Xa ^(x a , x fe ; h) , (2.12) 

6=1 \ n V X b)) 



and 



[V 2 /(x a )] = £ #iv^(x a; x 6 ; (2.13) 

respectively. The SPH equations, in their sum representation, can be directly derived by 
substitution of ( j2.1Q|) -( 2.12|) into the hydrodynamic conservation laws written in a Lagrangian 



form (Benz 1990). 

Higher order derivatives can be obtained following a procedure similar to that used to 
derive Eqs. ( |2.10|) - (|2.11|) . In principle, the only requirement needed is that a kernel of class C n 
is used in the smoothing of a n-order derivative. For instance, the smoothed approximation 
to the Laplacian operator, given by Eq. ( |2.13| ) requires a C 2 kernel, so the popular spline 
kernel (Monaghan & Lattanzio 1985) used in SPH fails because of its discontinuous second 
derivatives. For this reason, we will use a gaussian kernel 

where u = |x a — Xf,|//i, d = 1, 2, 3 denotes the dimension of the problem, and u a the kernel's 
compact support radius (Figure 1). As we will see below, a disadvantage of this kernel is 
that it needs to be defined over a larger compact support, typically u Q ~ 5. In contrast, the 
spline's compact support radius is by definition u Q = 2; hence, in general, using a gaussian 
kernel requires a larger number of neighbors to achieve similar accuracy as with the spline 
kernel. 

It is important to notice that we have distinguished the smoothing procedure ( ) from 
the discrete smoothing procedure [ ] to emphasize that there are in general two sources 
of truncation errors in approximating a function for its smoothed-discrete value; one error, 
the smoothing-error ((e)), arises from the smoothing procedure itself (see Eq. ( |2.2|) ), and a 
second error, the sum-error ([e]), is generated from the discretization of the integrals; that 
is 

N firA d n d n 
g— J n{Xb) dx n dx 11 

It is straightforward to show that for the one-dimensional gaussian kernel, the generalization 
of the smoothing-error ( p.2|) to n-order derivatives is 

(e) = T^ /(x) + 0(/l3) - (216) 

The distribution of the particles plays a fundamental role in computing the truncation 
sum-error. If the particles have a chaotic distribution, the approximations ( |2.10| )-( [2.13| ) 
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would be Monte Carlo estimates of the integrals in Eqs. (|2.1|) and ( p.6| )-( p78|) , respectively. 
However, in SPH the particles do not follow a chaotic position distribution; they are not 
completely independent. On the other hand, they are not, in general, distributed in a 
regular configuration either. Nonetheless, for estimating sum-errors, we are going to assume 
a regular distribution of particles because it is in this case for which the minimum truncation 
sum-error is obtained. This will guarantee that features arising from the sum-error in regular 
particle distributions will also show up in more realistic configurations. For simplicity, we 
will only consider here the one-dimensional case. 

Following a procedure similar to that of deriving truncation errors in standard quadrature 
formulas, it is not difficult to show that the leading truncation sum-error in ( |2 . 1 5|) is given 
by 

/Ar\ 2 d n+2 

[e] = h- n (-f) f(x)^- 2 w(u) + 0((Ax/hrh- n ), (2.17) 

where Ax is the particle interseparation and w{u) = hW(u,h) only depends on u. Fur- 
thermore, h/Ax = Ni/ 2 /u , where u Q is the compact support radius and N\/2 denotes the 
number of neighbors of a particle at each of its "sides" (Figure 1). It is important to notice 
from (|2.17 ) that in the case of "pure" smoothing (n = 0), the sum-error only depends on 



hj Ax oc iVi/2, the number of neighbors. 

Convergence tests in SPH have been usually done by doubling the total number of parti- 
cles; for our one-dimensional case, this implies 

Ax (k) = 2~ k Ax n with k = 0, 1, ... and Ax Q the 
spacing of the coarsest mesh. Little mention has been given, though, to how the number of 
neighbors is handled, i.e., the hj Ax oc iVi/2 — > oo limit. One has to remember that there are 
two competing errors in going from the continuum expressions to the discrete, smoothed ap- 
proximations. The smoothing-error has a discretization scale given by the smoothing length 
h, see Eq. fl2.1Gp , and the sum-error depends, on the other hand, not only on the number of 
neighbors hj Ax oc N\/2 within the compact support of the kernel, but also on the smooth- 
ing length h for terms involving differentiation, see Eq. Q2.17 ) with n > 0. Convergence 



tests have to satisfy the condition h~ n (Ax/h) 2 — > as Ax and h — ► 0. To illustrate this, 
we have computed the smoothed values of an analytic function and its smoothed first and 
second derivatives for different resolutions; that is, Ax^ k ' = 2~ k Ax Q and hP> = 2~ l h Q with 
k, I — 1, ..,4. Figure 2 shows the results using f(x) = sin(x). The lines join points of con- 
stant Ax. As expected, the top two plots show the second order convergence rate of |/(x)]. 
A similar second-order behavior is obtained for [/'( x )l an d [/"(^l if ^1/2 > 8. However, 
for N1/2 < 8, there seems to be an ill defined limit for [/'( x )l an d [/"(^)] when h — > 0; there 
the errors grow as h — > 0. As mentioned before, the answer to this puzzle is in how the limit 
h- n (Ax/hf is taken. We see that in particular the limit used in Fig. 2, Ax, h — > with 
hj Ax oc N1/2 constant, will break down for < 8. Of course in practice, it is impossible 
to take the limit to the continuum; however, it is important to have, at a given resolution 
scale, enough number of neighbors, so the smoothing procedure falls within the second-order 
regime of the truncation errors. In our one-dimensional case, this means to have N1/2 > 8. 

Finally, in SPH the particles are identified with fluid elements; hence, their dynamics is 
governed by the hydrodynamic equations. In generalizing SPH methods beyond the scope 
of hydrodynamics, we view SPH as strictly an interpolation scheme (Smoothed Particle In- 
terpolation: SPI) and, as such, used only to compute spatial gradients without the need of 
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introducing a grid. The particles could be either "fixed" in space, although not necessar- 
ily uniformly distributed, or free to move. The latter case will represent a truly adaptive 
implementation. For the moving particles (Lagrangian) case, a rule that dictates motion 
of the particles should be provided. In principle, the rule for moving particles can be arbi- 
trarily chosen; however, it is not clear at this stage whether non-single valued velocity fields 
are recommended since they could possibly lead to difficulties such as the one arising in 
hydrodynamic calculations when particle penetration occurs (Monaghan 1989). 

3. EULERIAN SPI 

In this section we will concentrate our attention on testing only the interpolating features 
of SPI. Hence, the interpolating particles will remain fixed in their original positions through 
out the evolution; the adaptive (moving particles) properties of SPI will be addressed in the 
next section. For these fixed particle tests, we start by considering the wave equation 

«9 2 0(x) - c 2 (x)V 2 0(x) = 0, (3.1) 

where d\ = J^-. We define the "momentum" variable x = dt<t> , and rewrite Eq. ( |3.1| ) as 

d t cj> = X . (3.2) 

and 

d tX = c 2 V 2 (3.3) 



We now apply smoothing to ( |3.2|) and (|3.3|) . It is important to keep in mind that the 
smoothing properties act only on spatial directions; hence, temporal differentiation and 
smoothing commute, (d t f) = d t (f). For temporal discretization, we use standard finite 
difference schemes. From substituting ( |2.10| ) and ( |2.13| ) into Eqs. ( |3~2| ) and ( |3.3| ), one obtains 
the discrete SPI approximations to Eqs. (|3.2|) and (|3.3|) to be 



.n+1/2 



At 

and 



xT y \ (3-4) 



n+3/2 _ v n+l/2 N in+l 

tt^ = c 2 a £ ^V 2 a W abj (3.5) 

At n b 

where W ab = W(x a , x b ; h). The time step is given by At = Dh/max(c s (x)) with D < y/3/2. 
It is evident from these expressions that temporal integration is performed using a leap frog 
method (Roache 1985). There are two interlaced time-steps (integer- and half-steps). The 
terms d t x, 4> an d spatial derivatives of (ft are placed at integer-steps; similarly, d t 4>, x an d 
spatial derivatives x are placed at half time-steps. This temporal discretization makes the 
numerical integration second-order accurate. 

To investigate the effects of the numerical dispersion in SPI, we evolved a gaussian pulse 
through a homogeneous medium; c(x) = 1 in Eq. ( |3.2| ). The initial conditions for and its 
"momentum" x were 

<p(x) = <p exp{-k 2 (x -t) 2 ] (3.6) 
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and 



X (x) = 2k 2 (x-t)4 > (x), (3.7) 
respectively. Solution errors are presented as cumulative relative errors; that is, 



(3i 



where and are the numerical and analytical solutions, respectively. Figure (3. a) shows the 
solution error as a function of the resolution scale following a displacement of four gaussian 
width-lengths. The second order convergence rate implied by Eq. ( |2.2| ) is evident from this 
figure. For comparison, in Fig. (3.b) we plot the difference between SPI and FD solution 
errors; the error differences are plotted as a function of resolution scale. The resolution 
scale is defined as the discretization scale (particle interseparation in SPI or mesh spacing 
in FD) for which both methods, SPI and FD, yield the same absolute truncation error. 
It is not difficult to show by substituting (|2.14j) into Eq. ( |2.2|) , that for uniformly spaced 



particles given a grid spacing Ax(FD) in FD, a smoothing length of h = Ax{FD)/\f?> yields 
the same SPI and FD truncation errors. Furthermore, one also gets that the ratio of SPI 
particle spacing to FD grid spacing is given by Ax(FD)/Ax{SPI) = y/3N 1/2 /u , where 
N1/2 is the number of neighbors at either side of each particle and u Q the radius of the 
support of the kernel in smoothing length units. As we mentioned before, typically u a ~ 5 
and N1/2 ~ 8; therefore, Ax(FD)/Ax(SPI) ~ 2.8. That is, for each FD grid point, we 
need almost three SPI particles to achieve the same absolute solution error. However, it is 
important to emphasize that both methods have the same convergence rate, namely 0(h 2 ) 
for SPI and 0(Ax 2 ) for FD. 

The next test constists of a one- dimensional gaussian pulse propagating through a medium 
in which the sound speed has a jump at the origin. That is, 

c(x) = c [3 + tan 2 h(x/(T)] (3.9) 

where, for simplicity, we use c Q = 1, and a is chosen such that the jump is smooth over a 
few particles. The pulse has a gaussian profile as before, see Eq. (|3.6|) , and is initially in 
the region x < where the sound speed is unity. Figure (4) shows four snapshots of the 
evolution. As expected, the pulse scatters at the point where the speed of sound has a jump. 
Two pulses emerge, one propagating in the same direction as the original pulse and a second 
one in the opposite direction. The calculation was repeated using a FD method; the bottom 
plot shows the difference between the SPI and FD solution for the last snapshot. 



4. LAGRANGIAN SPI 

We allow now for particle motions. However, for simplicity, we restrict ourselves to motion 
of particles in a constant velocity field. We start by rewriting equations Q3.2| ) and ( |3.3| ) in 
the frame of reference (t, x) comoving with the particles. The coordinate transformations 
between the comoving and laboratory frame (f, x') are t = t' and x = x' — vt', with v the 
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particle velocity. To avoid complicated notation, we have redefined the prime coordinate 
system as the laboratory frame. In these new coordinates, Eqs. (|3.2| ) and (|373|) take the form 



dtX ~ v ■ Vx = c 2 V 2 (4.1) 

and 

d t <j> - v • V0 = x, (4.2) 

respectively. One immediately sees that a description from a frame of reference comoving 
with the particles introduces transport-like terms in the equations; that is, terms of the 
form v ■ V but with the "wrong" transport velocity sign. These transport-like terms require 
the same careful consideration to prevent the development of instabilities as that given to 
transport terms in the density, energy and momentum equations in fluid dynamics. In 
addition to developing a stable form of the equations, a "correct" centering of each term 
in the equations must be implemented in order to preserve the second-order nature of the 
scheme. 

Let us look first at Eq. (|4.2|). As mentioned in the previous section, in a leap-frog 



integration method, spatial derivatives of V0 are placed at integer-steps; on the other hand, 
d t 4> and x are located at half-steps. Hence, an approximation to V0 at half time-steps is 
needed. We achieve this with the following second-order interpolation 

with similar interpolation for x; that is, it is necessary to save previous time-step values of 
cf) and x- 

Using (|2.10| ), (|2.11|) and (|2.13|) , one obtains the smoothed approximations of Eqs. (|4.1| ) 
and O to be 

j.n+1 _ xn N jn+l/2 N ,n+l/2 

• ' " ' 1 -v« • V a W ab + /3vlAtj: n V 2 a W ah + X a, (4.4) 



and 



At ^ n b ^ n b 



.n+3/2 _ n+1/2 N n+1 N in+l 



1^ = E —v a ■ V a W ab + CIJ2 —V 2 a W ab , (4.5) 



respectively. We have added to ( [4.1| ) a diffusion term, second term in the r.h.s. of Eq. 
of the form f3 v 2 At V 2 with (3 a dimensionless parameter of order unity. Without this 
numerical dissipation the equations are correctly centered but unstable. 

We repeated the simulation in Sec. 3 of propagating a pulse in a medium with a homo- 
geneous speed of sound, c(x) = 1 everywhere, but now the particle velocity is v = c; the 
pulse remains stationary in the particle frame. Because of the artificial dissipation term in 
Eq. (|4.4|) , the numerical solutions will drift from the analytic solutions of the continuum 
wave equation. Since we are interested for the moment in studying the convergence rate of 
SPI, we have added in Eq. ( |4.4| ) a source term —f3 v 2 At V 2 with <p the analytic solution. 
With this term, Eqs. ( POJ ) and ( ^4.5| ) have as analytic solution the gaussian pulse (|3.6|). This 
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allows a direct estimate of the solution error (Choptuik 1986) from Eq. (|3.8| ). We performed 
a series of runs with smoothing lengths = 2~ k h , k — 1, ..,4, where h a is the smoothing 
length of the coarsest level. The convergence factor of the solution error at the k refinement 
level is defined as (Choptuik 1991) 

£(fc) = _ (4 

where is the solution error (See Eq. |3.8| ) at level k. Figure (5. a) shows the convergence 
factor £( fc ) of the solution error as a function of time. The solution error obtained at the 
end of the runs after the pulse has traveled 3 width-lengths (t « 60) is plotted in Fig. (5.b). 
The solution errors suggest a convergence rate 0(h 139 ). This convergence rate implies that 
£( fc ) — > 2.62 as h — > 0, which is consistent with the behavior of in Fig. (5. a). 

Finally, we considered again the case in which the speed of sound has a jump given by 
(P-9|); the speed of the particles was set to v = 1, and the source term — [3v 2 AtV 2 of the 
previous test was not included. Figure (6) shows snapshots at the same time intervals as in 
Fig. (4); for comparison, in last plot of Fig. (6), solution differences between the moving and 
fixed particle cases are presented. 



5. SPI AND BURGERS EQUATION 

It is important to study the stability behavior of SPI discretization for nonlinear equa- 
tions. For this purpose, we consider now Burgers equation, which reads 

d t <P + ( t>d x <P = adl<p. (5.1) 

Due to the parabolic nature of Burgers equation, the temporal integration used in previ- 
ous sections to handle the wave equation becomes unstable when applied to this case. For the 
sake of simplicity, we will consider first the advection-diffusion equation which is obtained 
from Burgers equation after replacing <pd x <f) by vd x <p. Here again one needs to be careful 
with the correct centering of each of the terms in the equation. However, correct center- 
ing does not guarantee stability; it is well known that a centered-space and centered-time 
discretization of the advection-diffusion equation, although second-order accurate, is unsta- 
ble. A simple cure, at the expense of sacrificing second-order accuracy, is obtained by using 
forward-time, centered-space scheme. We are interested, nonetheless, in obtaining a scheme 
as close as possible to second-order accuracy and achieving, at the same time, stability. 

Our approach to solve the advective-diffusion equation is analogous to the FD scheme 
that combines leapfrog advection differences (0(At 2 , Ax 2 )) with forward-time, centered- 
space differencing of diffusion terms (0(At, Ax 2 )) (Roache 1985). For a oc O(At), this 
method is (0(At 2 , Ax 2 )). Under this method, the advection diffusion equation has a SPI 
approximation 

0a J a = -v* E fd Xa W ab + aJ2^-dlW ah . (5.2) 

As with the previous section, to study the convergence rate, we have added a source 
term to ( |5.2|) of the form —a d 2 ^, so an analytic solution to this modified advection-diffusion 
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equation is given by the gaussian pulse (|3.6|) . Figure (7. a) shows the convergence factor ^ 
of the solution error. In addition, the solution errors after the pulse has traversed three 
times its original thickness are plotted as a function of the particle separation in Fig. (7.b) 
for v = 1 and a = 1. These solution errors imply a 0(/i L95 ) scheme and — > 4 as h — > 0. 

A similar procedure to that used for the advection-diffusion equation was applied to 
Burgers equation. The stable SPI representation of Burgers equation is given by 

" J a = -<% E fd Xa W ab + aJ2^-dlW ab . (5.3) 

L b=l b b=l b 

Here again, to test the order of the method, we added a term d x (f) — a <9^0, so the equation 
has an analytic solution given by the pulse ( |3.6| ). The convergence factor of the solution 
error is plotted in Fig (8. a) and the final solution error in Fig. (8.b) for a = 0.2. We find that 
the SPI discretization (|5.3| ) to Burgers equation gives solution error of 0(/i 198 ). Finally, in 
Figure (9) we show the solution to Burgers equation, without the source terms (fid x (fi — ad x (j), 
for an initially gaussian pulse. As expected, the pulse develops shock-like features due to 
the non-linear term and also diffuses because of the dissipation term. 



6. CONCLUSIONS 

We have presented a method, smoothed particle interpolations (SPI), based on SPH 
techniques to solve hyperbolic and parabolic equations of non-hydrodynamic nature. The 
method was tested solving the one-dimensional wave equation with inhomogeneous sound 
speed and Burgers equation. We showed that the intrinsic adaptive nature of SPH is directly 
carried over. We also provided a prescription for handling advect ion-type terms that arise 
when particles are free to move or when dealing with Burgers equation. Furthermore, we 
showed the feasibility of obtaining 0(h 2 ) schemes under SPI discretizations. 

A possible drawback of this approach is its particle requirements compared with the 
mesh-point cost in finite difference methods. We have shown and tested that, for uniformly 
spaced particles, SPI required a minimum of three smoothing particles for each grid point 
in FD in order to achieve equivalent absolute solution errors. 

In principle, generalization of this work to two or three dimensional systems should 
be straightforward (Laguna 1994). After all, one of the reasons for SPH's popularity has 
been its simplicity for implementation in higher dimensions. Moreover, all the technology 
developed in N-body methods and standard SPH for finding neighboring particles can be 
directly utilized. 



ACKNOWLEDGMENTS 

I thank Richard Matzner and Matt Choptuik for numerous discussions and helpful sug- 
gestions. Work supported in part by the NASA (at Los Alamos National Laboratory), NSF 
Young Investigator award PHY-9357219, and NSF grant PHY-9309834. 



10 



REFERENCES 



Benz, W. 1990, in Numerical Modeling of Stellar Pulsation: Problems & Prospects, ed. J.R. 

Buchler, (Dordrecht: Kluwer Academic), p269 
Berger, M.J., & Oliger, J. 1984, J. Comput. Phys., 53, 484 
Choptuik, M. (1986) Ph.D. thesis, University of British Columbia 
Choptuik, M. 1991, Phys. Rev. D, 44, 3124 
Choptuik, M. 1993, Phys. Rev. Lett, 70, 9 

Davies, M.B., Benz, W. & Hills, J. 1991, Astrophys. J., 381, 449 

Davis, S., & Flaherty, J. 1982, SI AM J. Set. Statist. Comput., 3, 6, 

Evans, C. R. & Kochanek, C. S. 1989, Astrophys. J., 346, L13 

Hawley, J. F., Smarr, L. L., & Wilson, J. R. 1984, Astrophys. J., 277, 296 

Hernquist, L., & Katz, N. 1989, Astrophys. J., 70, 419 

Laguna, P., Miller, W.A., & Zurek, W. H. 1993, Astrophys. J., 404, 678 

Laguna, P., Miller, W.A., Zurek, W.H., & Davies, M.B. 1993, Astrophys. J., 410, L83 

Laguna, P. 1993, in preparation 

Mann, P.J. 1993, J. Comput. Phys., 107, 188. 

Monaghan, J. J. 1989, J. Comput. Phys., 82, 1 

Monaghan, J. J., & Lattanzio, J. C. 1985, Astron. Astrophys., 149, 135 
O'Rourke, P.J., & Amsden, A. A. 1986, Los Alamos National Laboratory report LA-10849- 
MS 

Roache, P.J. 1985 Computational Fluid Dynamics, (Albuquerque: Hermosa Publishers) 
Stellingwerf, R.F., & Peterkin, R.E. 1989, Mission Research Corporation report MRC/ABQ- 
R-1248 



11 



FIGURE CAPTIONS 



Fig.l Gaussian kernel with compact support radius u Q = 5. The compact support radius is 
defined as u Q = Ny 2 Ax/h where h is the smoothing length, Ax the particle intersepa- 
ration and A^/ 2 = 10 the number of neighbors of a particle at each of its sides. In the 
figure, the particle positions are denoted by open circles. 

Fig. 2 Smoothed values errors in the representation of f(x) = sin(x) and its smoothed first 
and second derivatives for Ax^ k ' = 2~ k Ax Ql k=l (circle), 2(cross), 3 (triangle), 4 (square) 
and h® = 2~ l h , I = 1,2,3,4. 

Fig. 3 Results from a gaussian pulse propagating in a medium with uniform sound speed 
c = 1. The particles were kept fixed throughout the evolution, (a) shows solution 
errors as a function of the particle interseparation. (b) shows the differences between 
SPI and FD solution errors. The errors were obtained after a pulse displacement of 
four initial gaussian width-lengths. 

Fig. 4 Propagation of a gaussian pulse through a medium in which the sound speed has a 
jump at the origin given by Eq. ( |3.9| ). The particles were kept fixed throughout the 
evolution. Bottom plot shows the difference between the SPI and FD solution in last 
snapshot. 

Fig. 5 Results from propagating a gaussian pulse in a medium with uniform sound speed 
c = 1. The particles moved with speed v = 1. (a) shows the convergence factor 
£( fe ) of the solution error, (b) shows the solution errors as a function of the particle 
interseparation at a pulse displacement of four pulse width-lengths. 

Fig. 6 Same as in Fig. (3) but for particles moving with speed v — 1. The plot at the bottom 
in this case shows the differences of the solutions between the fixed and moving particle 
cases at the end of the run. 

Fig. 7 (a) Convergence factor of the solution error for the advective-diffusion equation. 
A source term of the form —a d^cf) was added to the equation in order to have a 
gaussian pulse as analytic solution, (b) Solution errors as a function of the particle 
interseparation. 

Fig. 8 Same as in Fig. (6) but for Burgers equation with a source term d x (f) — a dl<f). 
Fig. 9 Solution to Burgers equation of an initial gaussian pulse with a = 0.2. 



12 



This figure "figl-l.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/9402062vl 



This figure "fig2-l.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/9402062vl 



This figure "figl-2.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/9402062vl 



This figure "fig2-2.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/9402062vl 



This figure "figl-3.png" is available in "png" format from: 



http://arXiv.org/ps/astro-ph/9402062vl 



